jst

preamble

library(scuttle)
library(BiocParallel)
library(SingleCellExperiment)
source("code/_utils.R")
bp <- MulticoreParam(3)
set.seed(251103)

loading

pbs <- readRDS("outs/qbs.rds")
sce <- readRDS("outs/fil.rds")
ist <- readRDS("outs/lv1.rds")
lapply(pbs, colnames)

$mye
[1] "DC"      "mono.c"  "mono.nc" "pDC"    

$nkt
[1] "NK/ILC" "Tc"     "Tcm"    "Tcn"    "Tex"    "Tfh"    "Thn"    "Tp"    
[9] "Treg"  

sub <- sort(unique(ist$clust))
(names(sub) <- sub)

[1] "mye" "nkt" "str" "tum"

selection

names(ks) <- ks <- unique(colnames(mx <- ist$profiles))
mx <- sapply(ks, \(k) rowMeans(mx[, k, drop=FALSE]))
sel <- lapply(ks[ks != "tum"], \(.) {
    my <- mx[, setdiff(ks, .)]
    fc <- rowMins(apply(my, 2, \(x) mx[, .]/x))
    fc <- fc[!is.na(fc) & is.finite(fc) & fc > 1]
    gs <- names(tail(sort(fc), 2e3))
})
sapply(sel, length)

str mye nkt 
841 729 360 

clustering

jst <- bplapply(sub[sub != "tum"], \(.) {
    ns <- c(1e4, 2e4, 1e5)/4
    nk <- seq(2, ifelse(. == "mye", 4, 8))
    if (is.null(gs <- sel[[.]])) gs <- TRUE
    se <- sce[, names(which(ist$clust == .))]
    ist <- .ist(se, nk=nk, ns=ns, gs=gs, pbs=pbs[[.]])
    ks <- intersect(letters, unique(ist$clust))
    ks <- c(ks, sort(colnames(pbs[[.]])))
    ist$clust <- factor(ist$clust, ks)
    ist
}, BPPARAM=bp)

kst <- ist
idx <- ist$clust == "tum"
kst$prob <- kst$prob[idx]
kst$logliks <- kst$logliks[idx, ]
kst$clust <- factor(kst$clust[idx])
jst$tum <- kst

visuals

lapply(setdiff(sub, "tum"), \(.) {
    kid <- droplevels(kid <- jst[[.]]$clust)
    sce$ist <- kid[match(colnames(sce), names(kid))]
    .plt_fq(sce, "sid", "ist", id=., h=TRUE)
}) |> wrap_plots(nrow=1)

ps <- lapply(setdiff(sub, "tum"), \(.) {
    kid <- droplevels(kid <- jst[[.]]$clust)
    sce$ist <- kid[match(colnames(sce), names(kid))]
    .plt_fq(sce, "ist", "sid", id=.)
}) 
ws <- sapply(ps, \(p) nlevels(p$data$Var1))
wrap_plots(ps, guides="collect", widths=ws)

mye

all

a

b

c

d

DC

mono.c

mono.nc

pDC

nkt

all

a

b

c

d

NK/ILC

Tcm

Tcn

Tex

Tfh

Thn

Tp

Treg

str

all

a

b

c

d

e

f

g

h

[[1]] NULL

[[2]] NULL

[[3]] NULL

labeling

lab <- list(
    str=list(
        epi="a",
        FRCcts="b",
        BEC="c",
        LEC="d",
        FRCtcz="e",
        FRCpv="f",
        fib="g",
        mCAF="h"),
    mye=list(
        tum="a",
        mac.pi="b",
        AML="d",
        DC="DC",
        pDC="pDC",
        TAM="c",
        mono.c="mono.c",
        mono.nc="mono.nc"),
    nkt=list(
        Tha="a",
        `NK/ILC`="b",
        Tp="Tp",
        Tex=c("Tex", "c"),
        Tfh="Tfh",
        Treg="Treg",
        Tcn="Tcn",
        Tcm="Tcm",
        Thn=c("Thn", "d"),
        NK="NK/ILC"))
kst <- lapply(names(lab), \(.) .jst(jst[[.]], lab[[.]]))
kst <- lapply(kst, \(.) {.$clust <- factor(.$clust, exclude="x"); .})
names(kst) <- names(lab); kst$tum <- jst$tum
sapply(kst, \(.) table(.$clust, exclude=NULL))

$str

   BEC    epi    fib FRCcts  FRCpv FRCtcz    LEC   mCAF 
  1034    212    653    941    744   3753    945   2016 

$mye

    AML      DC  mac.pi  mono.c mono.nc     pDC     TAM     tum 
   2943    1700    1529    1276     213    1587     866    1750 

$nkt

    NK NK/ILC    Tcm    Tcn    Tex    Tfh    Tha    Thn     Tp   Treg 
    42    422    197    747   1969    921    310   3513    358   1071 

$tum

  tum 
60121 

kid <- lapply(kst, \(.) .$clust)
idx <- names(kid <- unlist(unname(kid)))
kid <- kid[match(colnames(sce), idx)]
pbs <- aggregateAcrossCells(sce, kid[idx])
pbs <- logNormCounts(pbs)

n <- 7
ks <- colnames(pbs)
es <- logcounts(pbs)
gs <- lapply(ks, \(i) {
    j <- setdiff(ks, i)
    fc <- apply(es[, j], 2, \(.) es[, i]/.)
    fc <- rowMins(fc)
    names(tail(sort(fc), n))
}) |> unlist() 
mx <- .z(t(es[gs, ]))
pal <- setNames(.pal_kid[seq_along(ks)], ks)
pheatmap::pheatmap(mx, annotation_legend=FALSE,
    legend=FALSE, show_colnames=TRUE,
    treeheight_row=5, treeheight_col=5,
    fontsize=8, cellwidth=8, cellheight=8,
    cluster_rows=FALSE, cluster_cols=FALSE,
    color=pals::coolwarm(), border_color=NA,
    gaps_col=cumsum(rep(n, length(ks))),
    annotation_colors=list(x=pal),
    annotation_names_col=FALSE,
    annotation_col=data.frame(
        row.names=gs,
        x=rep(colnames(pbs), each=n)))

appendix

saving

saveRDS(jst, "outs/jst.rds")
saveRDS(kst, "outs/lv2.rds")

session

sessionInfo()

R version 4.5.1 (2025-06-13)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.6.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Madrid
tzcode source: internal

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] scran_1.36.0                patchwork_1.3.2            
 [3] ggplot2_4.0.0               tidyr_1.3.1                
 [5] dplyr_1.1.4                 BiocParallel_1.42.2        
 [7] scuttle_1.18.0              SingleCellExperiment_1.30.1
 [9] SummarizedExperiment_1.38.1 Biobase_2.68.0             
[11] GenomicRanges_1.60.0        GenomeInfoDb_1.44.3        
[13] IRanges_2.42.0              S4Vectors_0.46.0           
[15] BiocGenerics_0.54.1         generics_0.1.4             
[17] MatrixGenerics_1.20.0       matrixStats_1.5.0          

loaded via a namespace (and not attached):
 [1] gtable_0.3.6            xfun_0.54               htmlwidgets_1.6.4      
 [4] lattice_0.22-7          vctrs_0.6.5             tools_4.5.1            
 [7] parallel_4.5.1          tibble_3.3.0            cluster_2.1.8.1        
[10] BiocNeighbors_2.2.0     pkgconfig_2.0.3         pheatmap_1.0.13        
[13] Matrix_1.7-4            RColorBrewer_1.1-3      dqrng_0.4.1            
[16] S7_0.2.0                pals_1.10               lifecycle_1.0.4        
[19] GenomeInfoDbData_1.2.14 compiler_4.5.1          farver_2.1.2           
[22] statmod_1.5.1           bluster_1.18.0          mapproj_1.2.12         
[25] codetools_0.2-20        htmltools_0.5.8.1       maps_3.4.3             
[28] yaml_2.3.10             pillar_1.11.1           crayon_1.5.3           
[31] limma_3.64.3            DelayedArray_0.34.1     abind_1.4-8            
[34] metapod_1.16.0          locfit_1.5-9.12         rsvd_1.0.5             
[37] tidyselect_1.2.1        digest_0.6.37           BiocSingular_1.24.0    
[40] purrr_1.1.0             labeling_0.4.3          fastmap_1.2.0          
[43] grid_4.5.1              colorspace_2.1-2        cli_3.6.5              
[46] SparseArray_1.8.1       magrittr_2.0.4          S4Arrays_1.8.1         
[49] dichromat_2.0-0.1       edgeR_4.6.3             withr_3.0.2            
[52] UCSC.utils_1.4.0        scales_1.4.0            rmarkdown_2.30         
[55] XVector_0.48.0          httr_1.4.7              igraph_2.2.1           
[58] ScaledMatrix_1.16.0     beachmat_2.24.0         evaluate_1.0.5         
[61] knitr_1.50              irlba_2.3.5.1           rlang_1.1.6            
[64] Rcpp_1.1.0              glue_1.8.0              rstudioapi_0.17.1      
[67] jsonlite_2.0.0          R6_2.6.1